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ERROR PROBABILITIES FOR MAXIMUM LIKELIHOOD DETECTION 
OF M-ary POISSON PROCESSES IN POISSON NOISE 


By Sherman Karp and Martin G. Hurwitz 
Electronics Research Center 

and 

Robert M. Gagliardi 
University of Southern California 


SUMMARY 

In this problem we have considered some of the recent 
results in the detection of a Poisson-distributed signal in Poisson 
noise. Curves for error probabilities are presented for the case 
of detecting one of M equiprobable signals over a broad range of 
parameter values-. Implicit in these results for system applications 
is the use of "photon counting" receivers. Attention is given to the 
optical communication and radar problems for this receiver structure 
and significant parameters are translated into those used in the 
report. A complete description of the computational procedures 
used for making the error probability calculations is given. 


INTRODUCTION 


The purpose of this paper is to summarize some of the 
recent results (refs. 1-5) concerning the maximum likelihood 
detection of Poisson-distributed signals in the presence of 
Poisson-distributed noise and to tabulate the resultant error 
probabilities over a broad range of signal and noise when optimum 
signal design for maximum distance is used. Furthermore, the 
report shows how these results can be applied to the direct detection 
of optical signals, with the optimun • detector being a counter of 
photoelectrons. This form of detector can be implemented in the 
visible portion of the spectrum where photomultipliers exist. 

The tabulation is presented in two forms. The first form 
is related to the detection of M-ary signals and is applicable to the 
communications problem. The second form is related to the range 
bin problem in pulsed radar systems. It is felt that the values of the 
parameters will be applicable to most problems of this type. 



BACKGROUND 


When a classically describable field is incident upon a 
photodetecting surface, the probability dp of releasing a photoelectron 
in an interval dt over a surface area da is (ref 6): 


dp = al (t,a) dt da 

where I (t,a) is the intensity of the field and a is equal to Tj/hf in 
which is the quantum efficiency, his Planck's constant, and f is 
the frequency. For a surface of area A: 

O' dt / da I (t,a) = o'P (t) dt 
A 

where P (t) is the total collected power at time t. Suppose that the 
power is related to the incident particle rate n (t) by 

P (t) = n (t) hf. 

The probability of releasing a photoelectron in a time dt is then 
O' P (t) dt = Tjn (t) dt. 

With this assumption, the probability of releasing K photoelectrons 
p(K) in a finite interval AT is (ref. 7): 

[ r]n{t) dt]^ t + At 

p(K) = 1^7 exp [ - nn{t) dt] 


A system of events obeying this density distribution is called 
a Poisson process (ref. 8). If the system of events is stationary, 
we can replace the time average by an ensemble average, as 

p(K) = ^ exp[-? 7 nAT]. 


The system of events would then be called a stationary Poisson 
process (ref. 8). 

Strictly speaking, p (K) is a conditional density and should be 
written as p[K/n(t)] since, in general, n(t) itself is a sample from 


2 



a random process. However, in communications, one is interested 
in ’’designing" the waveform n(t) to satisfy certain desirable 
features. Consequently, n(t) is assumed to be deterministic. 

Let us assume that we are monitoring the current output of 
a unity quantum efficiency photodetector in an interval (0, T) and 
can distinguish all events. Let us further assume that one of two 
different rates was sent, resulting in one of two receive d n, (t), 
n^Ct), plus a stationary additive constant rate n corresponding to 
noise. We wish to formulate the maximum likelihood detection 
procedure for determining which rate is imbedded in the received 
signal. If the possible transmitted rates are both band-limited to 
B, then so are n^(t) and n 2 (t). We can therefore partition the (0,T) 
interval into M subintervals t-,, - t- = AT (ref. 1) 

1 J-T" i ^ 

where 


0 = t < t- < . 
o 1 


< t. 


M 




with AT consider the quantized version of the possible 

rates. That is; 


n. . 


1 

AT 


Li 


t. 

1 


n. 

J 


(t) dt; 


(ti 


<t ^ti^ 


). 


j = 1,2. 


To accomplish the_detection we consider two hypotheses: 
rate n^(t) and noise n^^ are present; H 2 “-rate n 2 (t) and noise 
are present. We now consider a vector space of M dimensions 
where each dimension represents the number of events observed 
in the corresponding interval. Since the number of events is 

independent from interval to interval, the vector K = (K-^ • . > 

K]^) has a conditional probability given of 

n )ATp 
” !— e 


M 


p(K IH.) = n 
i= 1 


[ (njj + 


K 
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and a .'onditional probability, given H 2 > of 


M 


K, 


p(K 


IH 2 ) = II _ - (n2i + nn) ^ 


i=l 


K. ! 
1 


The likelihood ratio A is then defined as 


p(kIh,) 

A(K) = ^ 

p(KiH2) 


M 

II 

i=l 


(n, . + n ) 

(n^. + n ) 

' Zi n J 



n2i)AT 


and the maximum likelihood detection criterion requires, after ob- 
serving K , a comparison of A(^) to a threshold c . If A^ c, we 
decide rate ni(t) is imbedded within the received process, and if 
A <c, we decide rate n 2 (t). Since the log function is monotonic, 
we can also make a decision based upon 

log A ^ log c 

where < denotes the threshold comparison. 

Thus log A is 


M 

i=l 


K. log 



+ n 


n 



n 

n 


M 

s - n^,^) A T ^ log c, 

i=l 


which can be rewritten as 


M 

A’ < [log c + S (n^^ - “2i» ^1 


( 1 ) 


where 
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( 2 ) 


M n. . M . 

A’ = 2 K, log (1 + - 2 K. log ( 1 + ) . 

i=l n i=l ^ n 

n n 

Thus the maximum likelihood detection te st inv^clves, equivalently, a 
comparison of the quantity A’ in Eq. (Z) to the new threshold in 
Eq. (1) . 

"Distance” Considerations 

Let us first consider the case where n^^ = 0 for all i and 
we are detecting the rate n^(t) + n^ versus n^ alone (i.e. , signal 
to no signal). Then the test in Eq. ( 1) becomes 

M n. . M 

A' == 2 K. log ( 1 + ~) ^ log c + 2 . ZTT 

i=l n ^ i= 1 

n 

Clearly, since is Poisson-distributed: 

M _ n-. . 

\ A'/h,] = S (n + n ) log (1 + — ) AT 

and 

M __ n, . 

Ej^[A'/H 2 ]= 2 log (1 + ^)AT . 


We define the "distance" D between the two hypotheses as 
D = Ej^[A'/Hj - [ A'/H^l 

M n 

- 2 n,. log (1 + ) AT . 

• 1 li 


We would like to investigate maximum and minimum of the distance 

M n . 

D = 2 log f 1 + zr-^] AT 

i=l 
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subject to the energy constraint 


M K 

S n = K' = — — . 

i=l AT 


Using a Lagrange multiplier we seek the minimum of 


M 

I = AT 2 n,. 
i=l 


n M 

log [ 1 + + X 2 

n i=l 

n 


n 


li 


which requires 

n,. AT n, . 

= AT log [ 1 + ^ + X = 0 

9^1i ^n ^i + ^n 

and is a set of M equations that must be satisfied for all i . 
solution is n^^^ = n for all i yielding 


M 


^li 


D . = 2 n, . log r 1 + ] at 

min . , li — 

1=1 n^ 


K 

= AT Mn log r 1 + — ] = K log [ 1 + —£ ] 

n ® Tn 


n 


n 


which can be identified as a minimum. 

The maximum value can be obtained by noting that 


M n.. M 2 

2 n log [1 + ^ [ 2 (n ) ] 

i=l n^ i=l 


M 

2 

i=l 


- n, . 

2 log^ (1 + z^) 


n 


1/2 


n 



The 


6 



with the equality occurring for = K ' . 

The maximum value for the sum becomes: 


D 

max 


M 


^li 


K 


S n log r 1 + ^ ] AT = K; log [ 1 + ] 


i=l 


n 


n 


n AT 
n 


Thus, the distance D is maximized when the rate is concen- 

trated in one AT interval. Hence, the signal maximizing detect- 
ability* (refs. 1,2) also maximizes distance. Notice also, that for 
the important case when for all i , 

M 2 
S n, . AT 


n 

n 


with distance and detectability identical. 


The maximization can also be obtained for the 


case n 21^0 . 


M . + n 

D = A'/hJ - E^, [ A'/H^] = S (n^. - log ( _^ ) AT 

i=l n^. + n 


2i 


n 


M 

i=i 


ni . + n 
li n 


, n^ . + n 
L 2 i n 


^^li ■ 


AT 


Since n^^ and n^^ must be non-negative, the log term is 
maximized for all i by having ri 2 j_ ~ 0 when nj^^ is not, and vice 
versa. This maximizes 


ni • + n 
^ 3 . n 

n^. + n I 
2 i n J 


(nii - n^.) 


by giving the largest magnitude to both the bracketed term and the 
exponent simultaneously and leaves the sign of the log positive. 

* Signal-to-noise ratio 
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Therefore, the two signals should be disjoint in time. On the 
assumption that they are disjoint, then; 


M n M 

D = S n, . log f 1 + 3 !^] + S 

i=l ‘ U1 


"2i f ^ 1 

n 

n 


To maximize D, we have only to notice that this is identical to 
maximizing each signal independently or to concentrate each signal 
in a different time interval. The optimum processor, therefore, 
calculates 


M lly M 

A’ = 2 K. log r 1 + 3 ^] - S K. log [ 1 + — 


n^. 

Zln 


i=l 


n i=l 


n 


n 


n 


which is compared to the threshold 


M 


log c + 2 (n.,. - n^. ) AT . 


i=l 


If the two signals are equiprobable with equal energy 


M M 

2 n, . AT = 2 n... AT , 

1=1 1=1 

c= land the threshold is zero. The H,, log likelihood functions 
are calculated for each choice of wavelorm n^ (t) and the hypothesis 
is selected according to the largest result. 

If M waveforms are used, the optimum processor would 
calculate the likelihood function for each of the M choices and 
select the maximum. Expressed mathematically, one obtains for 
optimum signal design under an equiprobable, equal energy assump- 
tion 


A' = K. log r 1 + 
1 


K 


n 


AT 


n 



log [ 1 + 


K 

<>o 

n AT 
n 
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or K. > K. . 

1 < J 

Hence, the problem is reduced to counting the number of photo- 
electrons in each AT interval and selecting the interval with the 
largest count. The probability of correct detection is then; 


= r Probability that K. > K. for all i/n. is the transmitted 
D ' r 1 J 1 J 

waveform]. 


M 
+ 2 
r=2 


r 


I" Probability that Kj = for r-l intervals/m is the 
transmitted waveform]. 


This is then averaged over all choices of waveform, which for 
the equiprobable set is just P^ • 

This can be written 


oo 

X=1 


(K + nn AT) _ + n AT) 

e s n 

(x) ! 


lM-1 


x-1 (n AT) .rr,\ 

^ ^ g(nnAT) 

i=o i ! 


M 

(1 + B) - 1 


MB 


,1 - (K + n M AT) 

+ — e ' s n ' 

M 


( 3 ) 


whe re 


B = 


(n_^AT) 


X 


X ! 


x-l (n^ AT)" 


1=0 1 i 


The error probability = 1 - 
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Presentation of the Data 


The data are presented in two forms. In the first form 
(Program 1, Figures L^IO) P^^ is plotted versus M for various 
values of Kg and Kn = AT. These parameters are directly 
related to the received signal and noise energies, respectively- 
Thus, if one determines from the range equation that Pg is the 
received signal peak power, then 


K 


s 


T] P AT 
' s 

hf 


whereas if the received average noise power is P^, 


K 


n 


hf 


(4) 


If non-diffraction-limited collecting optics is used, P is 
calculated as 


P., = A^T N. A\ 
N R o r X 


(5) 


with 


N. = spectral radiance (power/unit area, solid-angle band- 
width) 

T^ = fraction of optical transmission through all elements 

A = area of collector 
R 

= resolution of receiving collector (solid angle) 

AX= optical bandwidth. 

o 

In the visible region of the spectrum 5-6000 A, for non-dif fraction- 
limited optics with a blue sky background can be written as 

= 877AT(o:D)^ (AX) 
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_PROGRAMi 

_ 0^/30/68 80/80 LIST 

HURWl »C76501 »T500»CM50000. NASA I0l074 KARP3 HURWITZ 
REQUEST TAPE7.LO. PLOT SAO 568 RING IN. 

RFL. 50000. 

RUN(S'» .... .77777) 

LGO. ~ ~ ■ 


PROGRAM KARP3 ( INPUT. OUTPUT. TAPE7.TAPE5=INPUT.TAPE6=OUTPUT7 

DOUBLE PRECISION A , B . EKSN . ELI 0 . ELM . EM . EXK » PO . PE . P I ,PX . 

1 S. SAVE. SUM'. SUMP I . T . TERM . LRM2 , TERM3 , ftST . TOTAL. Xi . 

2, XK.XKN.XKNMAX.XKNMIN .XXS.XKSN.XM.XMl .XMB.XMBMIN.XX 

OATA IN. NOUT.IPLT. TEST/5, 6, Til. 0-24/ ' 

DATA XAX.YAX/8. ,11./ 

■ ■ DATA IXMAX,XKNMAX.XKNMIN,XMBMTN/li300.3. 0+2,1. D-30.1.D-12/ 

CCC KARP/HURWITZ 

CCC A'=ABS. VALUE OF OlFFERENCE BETWEEN TWO CONSECUTIVE I tRMS 

CCC B OESCRIBEO IN WRITE UP 

CCC ■■ EKSN = EXP (-XKSN) - - . . 

CCC EL10=CONSTANT USEO FOR TESTING MAGNITUDES OF NUMBERICAL VALUES 

CCC ELM = LOG(M) ' " 

CCC £M=INITIAL VALUE OF M 

CCC EX=ARRAY FOR LOG(M) FOR EACH KS AND EACH K 

CCC EXK=INPUT OF K 

CCC I = INDEX ■■■ ■ • 

CCC ICUTOF USED TO INSURE NO PREMATURE CUT-OFF OF THE CALCULATION 

CCC IEND = TEST VALUE FOR STOPPING RUN ' ' 

CCC lEND MUST BE NINE FOR ALL BUT LAST DATA CARD 

CCC ixmax=max value for index IX ■ ■■■ 

CCC 11 counts number total numbers of PE'S 

CCC J=lNDEX . . ... 

CCC jmp must be zero for first data card 

CCC “ jmp = tes't ‘ value for reading values of k 

CCC NM=NUMBER OF VALUES OF M 

CCC' NK = NuMBER OF value's O'f'k 

CCC NN=NPT(I) 

CCC N'PT (J) "counts 'number of M'S for' each k 
CCC N1=NN+1 

CCC ' 'N2 = NN + 2" ■ - 

CCC PD DESCRIBED IN WRITE UP 

CCC "'PE=l-'pb' 

CCC P I = ( KN**I) *EKSN/ ( I FACTORIAL) 

CCC PX = XKSN-»EKSN 

CCC S USED TO TEST MAGNITUDES 

CCC SAVE LAST VALUE OF TERM 

CCC SUM=SUMMATION (KN**I/I) 

CCC" S'UMP"l'= SUMMATION (PI) 

CCC T=DUMMY 

CCC TERM = PX*TERM2*TERM3 '' 

CCC TERM2=SUMPI**(M-1 ) 

CCC ' TERM3= ("l+B)**M-l/( M*B) 

CCC total=summat ION ( term ) 

CCC WY=ARRAY FOR LOG (PE) FOR EACH KS ANd’ EACH K 

CCC X=ARRAY FOR LOG(M) FOR ALL CURVES FOR EACH KS 

CCC XK=EXK(J) 

CCC XKN=KN 

CCC XKNMAX = MAX VALUE OF KN 

CCC XKNMIN=MIN VALUE OF KN 

CCC XKS=KS 
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80/80 LIST 


04/30/68 

CCC XKSN=KS+KN 
CCC XM=EM**I 
CCC XMBi=XM*B 
CCC XMBMIN=TEST value for XMB 
CCC XM1=XM-1 

CCC XX=IX=INDEX OF OUTER LOOP 

CCC Y=ARRAY FOR LOG ( PE ) FOR ALL CURVES FOR EACH KS 

DIMENSION X( 3000) >Y(3000) »EX( 100) »WY( 100) »NPT( IQQ) > EXK( lOO) 

EL10=300.DO*DLOG( lO.DO ) 

CALL INITPLTI IPLT) 

1 READ ( IN»3) 'EM»XKS. NM »NK ♦ I END ♦ JMP 
•- = 0 

IF (JMP .NE.O) GO TO 2 

READ (IN>4) (EXK(I)» I = 1 »NK) 

4 FORMAT (DIO. 3) 

CCC LOOP FOR K 

■‘2 DO 70 J=1*NK' 

NPTt J)=p 

3 FORMAT (2D10.3,4I5) 

XKN = EXK( J) 

CCC LOOP FOR M 
_ 10 DO 65 _I^ *NM_ _ 

CCC PRINT HEADINGS FTr'sT TIME THROUGH THIS LOOP 

IF ( I . NE.l ) GO TO 14 
WRITE(N0 u'T»12)XKS»XKN 

12 FORMAT ( 1H1«40X»3HKS=»D9,2>10X>3HKN=>D9.2// ) 

WRITE (N0UT.13) 

13 FO RMAT ( 17X>1HM,20X>2HPD>20X>2HPE>20X,6HLOG(M) » 16X_» 7HL0G ( PE ) / / ) 

14 XM=EM**f’ 

XM1=XM-1.D0 

CCC TEST TO insure EXP(KN) DOES NOT EXCEED MACHINE CAPACITY 
CCC IF KN GT XKNMAX» INCREASE M AND CONTINUE 

IF (XKN.GT.XKNMAX) GO TO 65 
XKSN=XKS+XKN 

CCC TEST VALUE TO INSURE NO PREMATURE CUTOFF 

ICUTOF=2.DO*XKSN 
EKSN=DEXP(-XKSN) 

CCC FORMATION OF PX FOR X=1 

PX=XKSN*EKSN 
SUM= 1.D0 

t=i.d6 

CCC FORMATION OF B FOR X=1 

B = XKN 

XMB = XM*B 

CCC TEST FOR (1+B)**M TOO LARGE 

S = DL0G ( 1 .D0+ 3) 

S=EL10/XM-S 

IF ( S.GT.O.DO) go to 17 

CCC value' for (1+B)**M too large for”machine 

TERM3=0.D0 

GO TO 16 

CCC VALUE F OR ( H- B ) **M SUFFIENTLY SMALL 

17 TERM3= ( (l.DO+B)**XM-l.DO )/XMB 

C(;C TEST FOR KN NEAR ZERO 

15 IF (x'KN.GE.XK'NMIN) GO TO 16 

TERM2=1«D0 

TERM3=1.D0 
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0^/30/68 . „ . 80/80 LIST 

GO TO 20 

CCC FORMATION OF P I ,SUMPI , TERM2 FOR X=1 
16 PI=DEXP(-XKN) 

SUMPI=PI 

TERM2=PI**XM1 

CCC LOOP FOR SUMMATION OF TERM 
20 TERM=PX*TERM2*TERM3 
TOTAL=TERM 
DO 50 IX=2»IXMAX 
SAVE=TERM 

xx=Tx 

X1=XX-1.D0 

CCC FORMATION OF PX FOR X.GE.2 
PX=PX*XKSN/XX 

CCC TEST FOR KN NEAR ZERO " 

IF (XKN.LT.XKNMIN) GO TO A5 _ 

CCC FORMATION OF PI FOR X.GE.2 

PI=PI«XKN/X1 _ 

CCC FORMATION OF TERM2 FOR X.GE.2 ' " 

CCC FORMATION OF B FOR X.GE.2 
SUMP I = SUMP I +P I 

25 TERM2= (SUMPI )**XM1 

T=t*XKN/Xl ‘ 

A = SUM 

SUM=SUM+T ’ - - ... 

B=XKN*B*a/ (XX*SUM) 

CCC TEST FOR (i+B)**M TOO LARGE 
S=DLOG(1.DO+B) 

s=Elio/xm-s ■■■ ” ‘ 

_ IF (S.LT.O.DO) GO TO 50_ 

CCC' FORMATION OF TERM3' 

CCC TEST FOR XMB LE XMBMIN 

CCC IF MB LE XMBMIN APPROX TERM3 WITH FIRST TWO TERMS OF MB ONLY 
35 XMB=XM*B 

IF IXMB.Gf.XMBMIN) GO TO 4b ' ' 

TERM3= l.DO+( XM-l.DO ) *B/2 .DO 

GO TO A5 

AO TERM3=((1.D0+B)**XM-1.D0)/XMB 
CCC FORMATION OF TERM FOR X.GE.2' 

45 TERM=PX*TERM2*TERM3 

CCC SUMMATION OF TERMS 
TOTAL=TOTAL+TERM 

CCC TEST IF DIFFERE'n'CE OF TERMS FOR X = N AND X = N + 1 IS SUFF. SMA'LL 
A=DABS( TERM-SAVE) _ 

IF (A.LT.TEST^AND. IX.GT. ICUTOF) GO TO 55 
CCC TEST TO INSURE AGAINST PREMATURE CUTOFF 
50 CONTINUE 

CCC. COMPUTATION OF FINAL VALUE FOR A GIVEN M 
55 PD=TOTAL+DEXP(-(XKS+XM*XKN) )/XM 
PE=1.D0-PD 
L = L+1 

ELM=DL0G10()!^M) _ 

X(L)=SNGL(ELM) 

ELM=DLOG10(PE) 

Y(L)=SNGL(ELM) ‘ 

NPT( J)=NPT( J)+l 

WRITE (NOUT»60) XM ,PD »PE »X ( L ) . Y ( L ) 
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04/30/68 80/80 LIST 

60 FORMAT(9X»3(D16,9,6X) ,2<F16.9»6X) ) 

65 CONTINUE 
70 CONTINUE 

L1 = L+1 

L2=L+2 

X (Ll ) = 0« _ 

X(L2)=0. 

Y(L1) = 0. 

Y(L2)=0. 

CCC PLOTTING RO UTINES FOR CAL-COMP PLOTTER 
CALL PLOTIO. .-31. .-3) 

CALL P LOT( 0. .2. .-3 ) 

CALL scale (X.XAX.L.l) 

CALL SCALE (Y.YAX.L.l) 

CALL AXIS (0. .0. .6HLOGIM) ,-6,XAX.6. .X(L1) »X(L2) ) 

CALL AX IS (0. .0. 7HLOG(PE) .7.YAX.90. . Y ( Ll ) . Y ( L2 ) ) 

11 = 0 

PO 90 1=1. NK 

NN=NPT ( I ) 

DO 85 J=1.NN 
11=11+1 

EX ( J ) = X ( 1 1 ) 

85 WY( J) = Y( iT) 

N1=NN+1 

N2 = NN + 2 - - . - 

EX(N1)=X(L1 ) 

EX(N2)=X(L2) 

WY(N1)=Y(L1) 

WY(N2)=V(L2) 

CALL LINE (EX. WY.NN. 1.0.0) 

90 CONTINUE 

CALL S YMBOL(0.5.16.0..3.3HKS=»0..3) 

CALL NUMBER( 1 .5. 16.0..3.XKS.O..-1 ) 

CALL PLOT( 20. .0. .-3 ) _ 

110 CONTINUE 

IF( IEND.EQ.9) GO TO 1 
CALL FIN(IPLT) 

120 STOP 

END 
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Figure 6,- Error probabilities for n^AT - fixed 
as a function of M 


■ ■ III I II I ■! 













where A\ is in angstroms, a i^ the resolution in arc seconds, and 
D is the diameter of the collector in centimeters. For diffraction- 
limited optics, a is related to D by 

5 X 

q: = Z.44 X 10 — arc seconds; 


hence 


K =4. 76^7 ATT (A\)\^xl0^^ 
n o 

with X the wavelength in centimeters. 

In the second form (Program 2, Figures 11-20) we consider 
the sampling interval T to be a constant and examine as M 

becomes large, or since AT = T/M, as AT becomes small- The 
assumption is that the signal energy can always be concentrated in 
the AT interval. For this case we use the parameters Kg and 
K = n^^T . The noise in the AT interval is then K/M . Kg can be 
calculated as before while K is calculated as 


K 


P T 
n 

hf 


It can be shown analytically that as M oo , an asymptotic 
value for P„ is reached, where (ref. 5) 


lim P^ 
M-^ oo 


1 + K 

c 





o 

K 


_ 

- 

_ 


This asymptote is quite apparent in Figures 11-20 and is seen to vary 
with K as expected. Physically, this implies that if bandwidth and 
computation are expedient, one can always approach 

= (1 + Kg) e"^s 

independent of the noise background by using narrower intervals. 
This procedure also provides better range resolution for the radar 
case. Implicit in these calculations is the fact that n^ is a constant 
or that the optical filter bandwidth is generally quite large. 
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PROGRAM 2 


04/30/68 '80/80 LIST 

HURW2,C76501 *T500»CM50000. NASA I0l074 KARP4 
REQUEST TAPE7»L0. PLOT SAO 766 RIfG IN. 



0000 . 

>♦»»»» 


77777) 


HURWITZ 


' PROGRAM' KARP4 ( I NPUT .OUTPUT . TAPE7 » TAPE5= I NPUT . TAPE6=0UTPUT ) 
DIMENSION EX( 100 ) .EXK( 100 ) .NPT ( 100 ) »WY ( 100 ) »X ( 3000 ) »Y( 3000) 
DOUBLE PRECISION A , B . EKSN . EL lO . ELM » EM . EXK »PD »PE » P I , PX . 

I S. SAVE, SUM, SUMPI.T. TERM. TERM2.TERM3.TEST.T0TAL.X1. 

2"“' Xk'.XKN.'XKNMAX.XKNMrN.XKS.XKSN.XM.XMl.XMBiXMBMIN.XX 

DATA IN.NOUT, IPLT,TEST/5,6,7,1.D-24/ 

DATA' XAX'.YAX/ll. jl4./' 

DATA IXMAX,XKNMAX,XKNMIN,XMBMIN/1000,3.D+2,1.D-30,1.D-12/ 

etc ' KARP/HURWITZ 

CCC A=ABS. VALUE OF DIFFERENCE BETWEEN TWO CONSECUTIVE TERMS 
CCC" 'B DE'SCRIBED IN WRITE UP 
CCC EXSN=EXP t-XKSN) 

CCC ELib='CONSTANT USED FOR TESTING MAGNITUDES OF NUMBERICAL VALUES 
CCC ELM=LOG(M) 

'CCC ' EM=INITIAL VALUE OF M 
CCC EXK=INPUT OF K 

CCC EX=ARRAY FOR LO'G(M') FOR' EACH 'kS AND EACH K 
CCC I=INDEX 

CCC I'l COUNTS NUMBER TOTAL NUMBERS ' OF' PE ' S 

CCC ICUTOF USED TO INSURE NO PREMATURE CUT-OFF OF THE CALCULATION 
CCC ‘lEND=TESf VALUE FOR STOPPING RUN 


.CCC lEND M.UST B£. NINE FOR ALL BUT LAST DATA CARD 
CCC IXMAX=MAX VALUE FOR INDEX IX 

CCC J=INDEX 

' CCC ■ JMP = TEST VALUE 'FOR READING VALUES' OF K 
__CCC JMP MUST BE ZERO FOR FIRST DATA CARD 

CCC nk=numBer of values' Of 'k 

CCC NM=NUMBER OF VALUES OF M 

'CCC' NN=NPT(I) 

CCC NPT(J) COUNTS NUMBER OF M'S FOR EACH K 
■ CCC ■ N1='nn+1 

CCC N2=NN+2 

Ce'e ‘ PD DESCRIBED' I N' 'W'R I T E UP ' 

CCC PE=1-PD 

CCC P 1= ( KN»*I ) *EKSN/ (I FACTORIAL) 

CCC PX=XKSN*EKSN 
CCC S USED TO TEST MAGNITUDES 
CCC SAVE LAST VALUE OF TERM 
Ce'e"' SU'M = SU'MMATI0N (kN**I/I) 

CCC SUMPI=SUMMATION (PI) 

CCC t=dum'my 

CCC TERM=PX*TERM2*TERM3 

CCC TERM2=SUMPI**(M-1) 

CCC TERM3= t 1+B)**M-1/ (,M*B),_ 

CCC ' T0TAL=SUMMATI0N ( TERM) 

CCC WY=ARRAY for LOG(PE) FOR EACH KS AND EACH K 
CCC X=ARRAY FOR LOG(M) FOR ALL CURVES FOR EACH KS 
CCC XK=EXK(J) 

CCC XKN=KN 

_ CCC XKNMAX=MAX VALUE OF KN 
'CCC' ■ '><KNM‘IN=MIN 'vALUE of KN 
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CCC XKS=KS 
CCC XKSN=KS+KN 
CCC XM=EM**I 

CCC X_MB = XM*B _ _ 

CCC XMBMIN=TEST VALUE FOR XMB 
CCC XM1=XM-1 

CCC XX=rX=INDEX ofouter loop 

CCC Y=ARRAY FOR LOG ( PE ) FOR ALL CURVES FOR EACH K.S 

EL10=300,Db*DLOG( lO.DO ) 

CALL INITPLK IPLT) 

r READ ( IN,3)’ EM,XI<;S, NM.NK. TEND. JMP ' ' 

L=0 

IF" CDMP '.NE.O) 'GO TO 2 ' ' 

READ (IN. 4) (EXK(l). 1 = 1, NK) 

4 FORMAT (DIO. 3) 

CCC LOOP FOR K 

2“DCr70“J=r.NK ■ ■ - - 

NPT( J) =0 

3 TORMAT r2Dr0.3,4T5) ■ 

XK=EXK(J) _ 

CCC LOOP FOR 'M 

10 DO 65 1 = 1, NM 

CCC 'PRINT' HEA'DINGS 'FIRST' TIME THROUGH THIS LOOP 
IF ( I.NE.l) GO TO 14 
■ WR'IT'E "(N0UT,12) 'XKS.XK' 

12 FORMAT ( IHl ,40X ,3HKS=,D9.2 , lOX ,2HK= ,D9.2// ) _ 

WRITE' (M0UT.13) ■ ■ 

13 FORMAT ( 17X,1HM,20X , 2HPD , 2 OX , 2 HP E , 20 X , 6HL0G ( M ) , 16X , 7HL0G ( PE ) / / ) 

14 XM''=EM**r 
XM1=XM-1.D0 

XKN = XK/'XM 

CCC TEST TO INSURE EXP(KN) DOES NOT EXCEED MACHINE CAPACITY 

CCC IF'“KN'GT "XK'NM'AX, rNCREASE M AND CONTINUE 

IF ( XKN.GT_.XKNMAX ) GO TO 65 
XKSN'=XKS'+XKN 

CCC TEST VALUE TO INSURE NO PREMATURE CUTOFF 
■ i'CUTOF = 2.D'0>XKSN 

EKSN=pEXP(-XKSN) _ _ 

CCC ' FORMATION OF'PX "FO'r' X='1 ' 

PX=XKSN*EKSN 

SUM=1.D0 

T=1.D0 

CCC ■ FORMATION OF B FOR 'X=l ' 

B = XKN _ _ . _ 

x'mb=x'm*b 

CCC TEST FOR (1+B)**M TOO LARGE 
S=DLOG(1.DO+B) 

S=EL10/XM-S _ 

IF (S.Gt.O'.DO) GO TO 17 ' ' 

CCC VALUE_FOR (1 + B)**M TOO LARGE _F(3R MACHINE 

TER'M3=0.D0 
GO TO 16 

CCC" ''V7n.ljFT’5R"(l+B)**M SUFFIENTLY SMALL 
17 TERM3= ( ( 1 .DO+B) **XM-1.D0 ) /XMB 
CCC ' ■ TES'T FOR KN NEAR ZE'RO 

15 IF (XKN.GE.XKNMJN) GO TO 16 

T'ERM2 = r.(56 
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. TERM3=1.D'0 

GO TO 20 

CCC FORMATION OF P I ,SUMPI , TERM2 FOR X=1 
16 PI=DEXP (-XKN) 

~ - - SUMPI=PI 

TERM2=PI**XM1 

"CCC "LOOP FOR SUMMATION OF TERM 

20 TERM=PX*TERM2*TERM3 

TOTAL=tERM ' " " 

DO 50 IX=2,IXMAX 

save=tErm — ■ — 

_ XX=IX 

X1 = XX-T.'D0' 

CCC FORMATION OF PX FOR X.GE.2 
'PX=PX*XKSN/XX 

CCC TEST FOR KN NEAR ZERO 

. . ..jp- (xi^^.LT.XKNMIN) go to 45 

CCC FORMATION OF_ PI FOR X.GE.2 

■“ Pr=Pl*XKN/Xl 

CCC FORMATION OF TERM2 FOI^ X.GE.2 

crc“" Formation OF B fott xVge.2 

SUMPI=SUMPI+PI 

25"TERM2=(SUMPI )**XM1 

T=T*XKN/X1 
A = SUM 
SUM=SUM+T 

B^‘XKN*B*A/(XX*SUM) - - - 

CCC TEST FOR (1_+B)**M TOO LARGE 

— ■ "s=dLogti.oo+B) ~ ■ ■ • - - - - 

S=EL10/XM-S 

IF 7 S.LT. 0 .D 0 I GO TO“5o 
CCC FORMATION OF TERM3 

CCC TEST"F0R XMB L£ XMBMIN 

CCC IF MB LE XMBMIN APPROX TERM3 WITH FIRST TWO TERMS OF MB ONLY 
'as ‘‘XMB=XM*B 

IF (XMB. GT. XMBMIN) GO_ TO 40 

- 1VD0+(XM-1.D0)*B/2.D0 

GO TO 45 

4(D ■TERM3= ( ( r.DO+B) **XM-r.OO ) /XMB 
CCC FORMATION OF TERM FOR X.GE.2 
45 TERM=PX*TERM2*TERM3 
CCC SUMMATION OF TERMS 
T0TAL=T0T AL+TE“RM 

CCC TEST IF DIFFERENCE OF TERMS FOR X=N AND X=N+1 IS SUFF. SMALL 

A=DABS( TERM-SAVE) 

IF (A.LT.TEST.AND. IX.GT. ICUTOF) GO TO 55 
■ CCC TEST TO insure AGArNST PREMATuRE CUTOFF " 

50 CONTINUE 

CCC COMPUTATION OF FINAL VALUE FOR A GIVEN "M 
55 PD=TOTAL+DEXP(-(XKS+XM*XKN) )/xM 
PE=1.D0-PD 

L = L+1 _ . 

ELM=DL0G10(XM) 

X(L)=SNGL(ELM) 

‘ ■elm=dlogio'(P'£ ) 

Y(L)=SNGL(ELM) _ 

NPf ( J)'=NPT('J) +I 
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WRITE (NOUT»60) XM,PDtPE,X(L) »Y(L) 

60 FORMAT (9X »3(D16.9*6xr^2<El6. 9 »6X) ) 

65 CONTINUE "" 

70 CONTINUE 
L1=L+1 
L2 = L,+2 
X(L1)=0. 

X(L2)=0. 

Y(L1)=0# 

Y{L2)=0. 

CCC PLOTTING ROUTINES FOR CAL-COMP PLOTTER 
CALL PLOT( 0. *-31. »-3 ) 

CALL PLOT( 0. *2. *-3 ) 

CALL SCALE (X*XAX>L,1 ) 

CALL SCALE (Y,YAX»L,ir‘ 

CALL AXIS ( 0, ,0. »6HLOG(M) ,-6*XAX»0*»X(Ll),X(L2)) 
CALL Axis' {O'. ,0»7HLOG(PE) ,T»YAX»90.»Y(L1) »Y(L2) ) 
I 1 = 0_ 

DO 90 ‘ r=T,NK * 

NN=NPT ( I ) 

DO 85 J=1.NN 
I 1 = 11 + 1 

EX(J)=X(II) ' 

85 WY( J) = Y( ID 

■ Nl=NN+r " ■■ ■ 

N2=NN+2 
hX(Nl)=X(Ll) 

EX(N2)=X(L2) 

WY(N1)=Y(TL1) 

WY(N2)=Y(L2) 

■ ■ CA-LU LINE'rEX,WY,NM7T ,0,0) 

90 CONTINUE 

CALL SYMBOL( 0.5 , 16 .O'r.3 ,3HKS= ,0. ,3 ) 

CALL NUMBER ( 1 .5 , 16 .0» . 3 ,XKS,0. ,~1 ) 

CALL PLOT( 20.', 0. ,-3 ) 

110 CONTINUE 

IF( i'END.EQ.9) GO TO 1 
CALL FIN( IPLT) 

120 STOP 
END 
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COMPUTATIONS 


The evaluation of and P„ as functions of the parameters 
Kg j and M, was accomplishea by a set of computer programs 

which are presented below. 

In one case. Kg, and M are input to the program. In the 

second case is computed as Kn = K/M where Kg, K and M are 
input to the program. 

The programs were written in Fortran IV and run on the 
CDC 6400 Computer System at the Smithsonian Astrophysical 
Observatory in Cambridge, Massachusetts. 
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DESCRIPTION OF THE FLOW CHART 


The constants including tape definitions are set in data 
statements (Figures 21a, b, and c): 


B ox 

1: 

Input data and switch values. Set index J to 1. 

This index counts the number of K values. 

n 

B ox 

2; 

Set index I to 1 . This index counts the number of 
M values. Print headings for tables to be outputted. 

Box 

3: 

Set the M value for this loop. In one version of the 
program K is fixed. In the other version it is 
computed as a function of M. 

B ox 

4: 

Test magnitude of to prevent overflow or error. 

B ox 

5: 

Establish initial values of computed parameters. 

B ox 

6: 

Test magnitude of (1+B)^ to prevent overflow or 
error. 

B ox 

7: 

Compute value of term 3. 

B ox 

8; 

Test for K near zero, 
n 

B ox 

9: 

Compute initial values for PI, Sum PI, Term 2. 

B ox 

10: 

Compute initial values for term, total. 

B ox 

11: 

Set index IX to 2. This counts number of times 
through major computation loop. 

B ox 

12: 

Compute PX. 

B ox 

13: 

Test for K near zero, 
n 

B ox 

14: 

Compute PI, Sum PI, Term 2. 

B ox 

15: 

Compute new B. 

B ox 

16: 

Test magnitude of ( 1 + B . 

B ox 

17: 

Compute Term 3. 
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B ox 


Compute new term and add this total, 


Box 

Box 
Box 
Box 
Box 
Box 
B ox 
B ox 
Box 

B ox 
B ox 


18: 

19: Test for completion of major loop. If the difference 

between two successive values of term is sufficiently 
small, and the index has exceeded a predetermined 
cutoff value, the loop is completed. 

ZO: Increment IX by 1 . 

21: Test for sufficient number of times through the loop. 

22: is set to final value of total. is l-P,^ . 

D ED 

23: Print out table of M, their logarithms. 

24: Increment index I by 1 . 

25: Test for maximum value of M. 

26: Increment index J by 1. 

27: If J has not exceeded J max, recycle for next 

value of K (or K ) . 

' n' 

28 : Plot re suits . 

29: Test whether any more data sets are to be computed. 



.START 


INITIALIZE 

INPUT 

I—^J 


I— ►! 

PRINT HEADINGS 

2 

EM^— ► M 
K/M — ►KN 


<J^>XKNMWf> 

KS+ KN — ►KSN 

(KSN)e”^^^ — ♦ PX 

I ► SUM 

KN-^B 


(I +8)'^ < 10^°^ 
^ 


(d+Br-D/MB 


e-KN ► PX 

PI ► SUMP I 

pltM-')__^TERM 2 


rKN<XKNMINJ 

\ ^ 


(PX){TERM2)(TERM3 
— TERM 


I 

1 — ► TERM 2 

TERM ► TOTAL 


1 — ► TERM 3 


Figure 21a.- Program flow chart 
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